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SUMMARY 


A  computer  program,  SPINFIT,  is  designed  for  determining  the  linear  and 
nonlinear  spin  damping  coefficients  by  fitting  calculated  values  to  the  observed  wind 
tunnel  spin  data.  The  previous  method  of  fitting  a  closed  form  (exponential  form)  of  the 
differential  equation  for  spin  rate  only  considered  linear  spin  damping  coefficients. 
The  results  from  this  exponential  fitting  were  found  to  be  unsatisfactory  for  some  cases 
of  the  wind  tunnel  test  data.  This  result  suggested  that  a  nonlinear  spin  damping 
coefficient  term  should  be  included  in  the  equation  of  motion. 

The  SPINFIT  computer  program  uses  an  equation  that  has  a  second  order  term 
in  spin  rate  to  extract  the  nonlinearity  in  a  spin  damping  coefficient.  This  was 
demonstrated  by  fitting  the  wind  tunnel  spin-up  test  data  for  the  81 -mm  M819  RP 
smoke  projectile.  The  result  shows  that  the  SPINFIT  more  closely  fits  the  data  when 
compared  to  the  exponential  fitting.  This  indicates  that  the  spin-up  test  data  has 
considerable  amount  of  nonlinearity  in  spin  damping  coefficient.  Because  of  the 
additional  second  order  term,  the  residuals  and  the  probable  errors  are  reduced,  and 
the  estimated  values  of  initial  spin  rate  and  steady  state  spin  rate  are  closer  to  the 
observed  wind  tunnel  test  data  when  compared  to  the  results  from  the  exponential  fit. 

The  SPINFIT  program  is  not  limited  to  the  reduction  of  wind  tunnel  test  spin 
data.  It  can  also  be  used  for  free  flight  test  data  when  sufficient  data  points  are 
available  for  the  fitting  process. 
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INTRODUCTION 


Exponential  fitting  is  done  by  matching  the  solution  of  the  linear  equation  of 
motion  to  the  observed  data  points.  This  method  assumes  that  there  is  no  nonlinearity 
in  the  spin  damping  coefficient.  Therefore,  this  method  may  be  used  successfully  only 
if  the  spin  data  contains  negligible  nonlinearities  in  spin  damping  coefficient. 
However,  in  many  cases,  the  large  nonlinearity  does  exist  and  the  second  order  term 
in  spin  damping  coefficient  must  be  added  to  the  equation  of  motion. 

A  computer  program  called  SPINFIT  has  been  developed  to  extract  the 
nonlinear  spin  damping  coefficient  by  adding  a  second  order  spin  damping  term  in  the 
differential  equation  of  motion.  Because  of  this  additional  coefficient,  it’s  able  to  match 
the  obsen/ed  data  curve  more  closely.  Even  though  the  first  order  differential  equation 
with  nonlinear  spin  damping  term  can  be  integrated  to  obtain  a  closed  form  solution, 
the  approach  used  in  this  program  can  handle  any  additional  higher  order  spin 
damping  terms  in  the  differential  equation  with  simple  modification. 

Chapman  and  Kirk  of  Ames,  National  Aeronautics  and  Space  Administration, 
developed  a  technique  which  is  used  to  extract  aerodynamic  coefficients  from 
projectile  angular  motion  by  numerically  integrating  the  second  order  differential 
equations  of  motion  and  fitting  the  results  to  the  observed  data  (ref  1).  This  report 
employed  this  technique  to  determine  the  linear  and  nonlinear  spin  damping 
coefficients  from  the  wind  tunnel  experimental  spin  data. 

ANALYSIS 

The  equation  of  motion  for  the  spin  of  a  test  model  is  assumed  to  be  in  the  form 
of 


qsd  f 

'■71 
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where  -L  is  bearing  friction  term  in  wind  tunnel  test. 
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So  rewrite  the  equation  in  terms  of  C’s 

P^c^+c^p+c^p  ^2) 

The  constant  C,,  C,,  and  can  be  determined  by  fitting  the  solution  of  the  assumed 
differential  equation  of  motion  directly  to  the  observed  wind  tunnel  spin  data  (P,). 

During  the  fitting  process,  the  initial  spinC,  and  the  spin  parameters 
Cj,  C3,  and  will  be  adjusted  until  the  best  fit  is  obtained.  The  best  fit  is  defined  as 
when  the  sum  of  the  squares  of  the  residuals  (Pj  defined  below  is  a  minimum. 

2 

/-I 

Where  N  is  the  total  number  of  data  points  and  t,-  is  the  time  at  the  ^  data 
point.  The  necessary  conditions  for  residual  to  be  minimum  with  respect  to  the  initial 
spin  C,  and  the  spin  parameters  C^,  C^.  and  C,  are 


SR. 


sc. 


^  BPAt,) 

=  (»,)-  /’c  (»/)l - = 00  k=1 ,2,3,4 


ac* 


(4) 


Since  P  is  an  implicit  function  of  Cy’s,  the  equation  has  to  be  solved  by  iteration 
procedure.  To  obtain  the  corrections  for  the  parameters  from  the  iteration  to  the 
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j+^st  iteration,  P  is  expanded  in  a  Taylor  Series.  Neglecting  the  second  and  higher 
order  terms,  then  the  result  becomes 

dPit,)  8 Pit.)  dP(t,)  BP  it.) 

l;^^4  (5) 

aCj  aCj  3C3  BC^ 

For  i=1 ,2 . N 

The  equation  (5)  is  substituted  into  equation  (4)  for  calculating  the  y+lsf  corrections. 
For  simplicity,  the  subscripts  j  and  j  -t-  1  will  be  omitted  in  the  following  expressions 
without  ambiguity. 


That  is 


2 


- AC, 

dC, 


where  k=1 ,2,3,4 


Since  the  summation  operator  is  reversible,  and  after  manipulation,  the  above 
equation  yields 


where  k=1 ,2,3,4 
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Rewrite  the  equation  in  matrix  notation 


[A1[C1  =  [B1 


-2( 

i-1 


l.k=1,2,3.4 


ac, 


ac. 


i-1  * 

Ck-^Ck  where  k=1 


k=1 .2,3,4 
,2.3.4 


Where  the  dimensions  are  (4  by  4),  (4  by  1),  and  (4  by  1),  respectively. 

Then  the  correction  of  spin  parameters  ACu  AC2,  AC3,  and  AC4,  can  be 
determined  by  [C  ]  =  [A  ]-i  [B  ]  if  and  only  if  the  [A  ]-i  exists. 

The  iteration  process  is  employed  to  obtain  corrections  of  the  spin  parameters 
from  the  equation  by  fitting  the  wind  tunnel  spin  data.  It  will  be  terminated  if  any 
prescribed  degree  of  convergence  or  number  of  iterations  is  reached.  The  method  of 
obtaining  the  partial  derivatives  is  to  integrate  the  parametric  spin  function  of  the  spin 
parameters  with  time. 


If  we  assumed  the  order  of  differentiation  can  be  reversed  then  that  is 


ac 


- }=- - 

ac,  dt  ^  dt'-  ac,  ' 


Then 


ai>(r,) 

Where  Qj(t,)= -  J  =  1,2,3, 4 

dCj 

G,=(C3+2c,p)ej 
Q2=(Cj+2C^P)Q^+1 
Q,=  iC,  +  2C^P}Q,+  P 
Q,=  iC,+2C,P)Q^+P^ 


with  (/,  )  omitted  for  the  P  and  Q’s 
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In  the  program,  equations  for  spin  and  Qj  calculations  are  integrated 
numerically  with  the  initial  conditions  of 


/’(0)=Ci 

Gi(0)  =  1.0 


G2(o)=G3(0)=e4(0)=o.o 


At  each  iteration,  the  corrections  for  the  spin  parameters  is  added  to  values  to 
get  new  spin  parameters,  and  the  new  spin  data  is  calculated  and  compares  with  the 
wind  tunnel  spin  data.  This  process  is  repeated  until  the  residual  reaches  its  minimum 
value  or  reaches  the  iteration  number  limit. 

RESULTS  AND  DISCUSSION 

The  SPINFIT  program  gives  improved  fitting  when  compared  to  the  exponential 
fit.  The  differential  equation  of  motion  for  spin  used  in  the  program  is  as  follows: 


isd  (  (Pd^  (Pd)-] 

-  C  8+C,  —  +C,  —  U- 

/  V  •  2V  r^\2V  ^  '  I 


Unlike  the  equation  above,  the  exponential  fit  results  from  the  solution  of  a 

differential  equation  which  does  not  contain  the  second  order  term  ^^p2^  The 
exponential  fit  cannot  take  into  account  the  nonlinearity  in  the  observed  data  because 
the  equation  doesn’t  have  higher  order  terms.  As  we  all  know,  more  higher  order 
terms  give  better  estimated  results,  but  in  most  cases,  terms  higher  than  the  second 
order,  are  not  needed  because  the  higher  order  terms  are  negligible  in  spin  damping 
coefficient.  The  SPINFIT  program  contains  a  second  order  term  to  handle  nonlinearity 
in  the  observed  spin  data  for  better  fitting. 

For  the  example  of  M819  RP  smoke  cartridge,  configuration  119123  (round 
137)  wind  tunnel  spin-up  test  data,  it  shows  larger  residual  from  the  exponential  fit 
than  from  the  SPINFIT  program.  According  to  the  computer  results,  the  spin-up  test 
data  has  nonlinearity  which  the  SPINFIT  program  is  able  to  detect.  The  calculated 
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spin  data  curve  from  the  SPINFIT  follows  the  test  data  more  closely  than  the 
exponential  one. 

The  sigificant  improvement  of  using  nonlinear  fit  method  to  calculate  initial  spin 
rate  and  the  steady  state  spin  rate  was  observed  from  the  example.  From  the 
computed  result,  SPINFIT  gave  0.05  rev/s  of  probable  error  for  the  initial  spin  rate 
where  the  exponential  fit  gave  0.12  rev/s.  Overall  residual  of  SPINFIT  result  was  0.33 
{rev/s) 2  when  the  exponential  fit  was  2.98  (rev/s)  2  .  The  SPINFIT  reduced  the 
probable  error  of  fitting  process  to  0.0098  rev/s  from  0.085  rev/s  of  the  exponential  fit. 
The  results  of  the  fits  are  summarized  in  tables  1  and  2. 

CONCLUSIONS 

The  nonlinear  and  exponential  fits  to  the  observed  data  for  81 -mm  M819  RP 
smoke  cartridge  spin-up  wind  tunnel  test  are  shown  in  the  figure.  The  fit  obtained 
using  the  nonlinear  spin  damping  term  in  the  equation  of  motion  is  improved  when 
compared  to  the  exponential  fit.  That  is  because  the  wind  tunnel  data  does  contain 
nonlinearity  and  it  is  handled  by  the  second  order  term  of  the  differential  equation  of 
motion  in  the  SPINFIT  program.  This  program  gives  excellent  estimations  of  an  initial 
spin  rate,  a  steady  state  spin  rate,  and  spin  damping  coefficients 
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Table  1 .  Linear  and  nonlinear  fitting  comparisons 


COMPARISON  j 

ITEMS 

EXPONENTIAL  FIT  }  SPINFIT  } 

I  I 

STEADY  STATE  SPIN  RATE  (rev/sec) 

29.331  i  26.7534  I 

1  1 

1 

SPIN  DAMPIN  MOMENT  COEFF  AT  STEADY  STATE 

-0.01096  j  -0.00866 

TOTAL  RESIDUAL  {rev /sec)' 

2.981  I  0.330  1 

1  i 

l_  _  _  -  _  . 

TOTAL  PROBABLE  ERROR  (rev/sec) 

0.2483  !  0.0846  i 

j 

1.34  j  1.89  j 

PROBABLE  ERROR  OF 

0.12  j  0.05  j 

1  1 

0.0008  1  0.0006  1 

1  1 

PROBABLE  ERROR  OF 

0.000015  j  0.000013  j 

c, 

P 

-0.011  1  0.004  j 

1  1 
. l._  .  1 

PROBABLE  ERROR  OF  C, 

c, 

0.000  0.000801 

1  1 

- ^ - 1 

N/A  j  -0.183  j 

1 

PROBABLE  ERROR  OF  C, 

N/A  1  0.010131  1 
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Table  2.  Observed  and  calculated  data  points 


OBSERVED  DATA  POINTS  AND  FITTED  RESULT  POINTS 


NUM. 

TIME 

OBSERVED  SPIN 

EXPONENTIAL  nT 

SPINFIT 

sec 

rev/sec 

rev/tec 

rev/sec 

1 

1.94 

1.82 

1.3363 

1.8898 

2 

2.41 

2.55 

2.3494 

2.6902 

3 

2.76 

3.24 

3.0799 

3.2881 

4 

3.05 

3.87 

3.6701 

3.7842 

5 

3.29 

4.44 

4.1486 

4.1949 

6 

3.51 

4.62 

4.5793 

4.5711 

7 

3.72 

5.0 

4.9836 

4.9299 

8 

3.92 

5.0 

5.3625 

5.2711 

9 

4.19 

5.85 

5.8647 

5.7307 

10 

4.51 

6.49 

6.4462 

6.2734 

“  i 

4.82 

6.67 

6.9959 

6,7965 

12 

5.10 

7.27 

7.4810 

7.2664 

13"! 

5.43 

7.83 

8.0392 

7.8163 

14 

5.80 

8.37 

8.6481 

8.4274 

15 

6.15 

9.00 

9.2081 

8.9993 

16 

6.48 

9.47 

9.7221 

9.5324 

;  17 

6.78 

10.00 

10.1781 

10.0114 

18 

7.08 

10.59 

10.6234 

10.4846 

19 

7.40 

10.91 

11.0870 

10.9827 

20 

7.76 

11.43 

11.5949 

11.5343 

21 

15.0  i 

20.25  ! 

19.2783 

20.1429 

22 

21.5 

23.90 

23.2928 

24.0043 

23 

25.5  1 

25.25 

24.9185 

25.1993 

24 

30.5 

26.10 

26.3497 

26.0038 

25 

32.5  j 

26.10 

26.7824 

26.1951 
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SPIN  RATE  (REV/SEC) 
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LIST  OF  SYMBOLS 


S 


d 

Lf 

lx 

Pe 

Pc 


q 


Qi 


s 

Cl 


C2 


C3 


Roll  moment  coeffient  due  to  fin  cant  (per  rad) 

Linear  spin  damping  coefficient  (per  rad) 

Nonlinear  spin  damping  coefficient  (  per  rad  2  ) 

Fin  canted  angle  ( rad ) 

Projectile  reference  diameter  ( ft ) 

Bearing  frictional  torque  ( Ib-ft ) 

Projectile  axial  moment  of  inertia  (Ibm-^^ ) 

Experimental  spin  data  ( rev/s ) 

Calculated  spin  data  ( rev/s ) 

Dynamic  pressure  ( \b/ft  2  ) 

Partial  derivative  of  spin  with  respect  to  Ct  ,62,03,  and  C4 
Reference  area  ( ff2) 

(  Pq  )  Initial  spin 


[  qsdCi,5 

I"  'x 


Constant  of  linear  spin  moment  term 


I  qsdCip,d' 


Constant  of  nonlinear  spin  damping  moment  term 


C4 


I  qsdCipd' 
[  2VI^ 


Constant  of  linear  spin  damping  moment  term 
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Vetocity  ( ft/s ) 

Free  stream  air  density  (slug/ft^) 
Derivative  ( )  respect  to  time 
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APPENDIX  A 

PROGRAM  SYMBOL  CONVERSION  TABLE 
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COMPUTER  CODE 

SYMBOL 

DIA 

d 

AX 

AXX 

DEL 

5 

FT 

AC(1) 

Po 

AC(2) 

C, 

i 

AC(3) 

c, 

AC(4) 

% 

QQ 

q 

PI(i) 

Wind  Tunnel  Spin  Data 

P(i) 

Computed  Spin  Data 

Q(N,M) 

e^(«) 

DT 

A  t 

TIME(i) 

Time  Table 

17 


APPENDIX  B 

COMPUTER  INPUT  FORMAT 
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CARD 


NUMBER 

COLUMN 

NAME 

FORMAT 

DESCRIPTIONS 

UNITS 

1 

1-72 

DIDENT 

18A4 

Identification  card 

- 

2 

1-5 

NT 

15 

Number  of  data  point 

- 

6-10 

DT 

F5.0 

Time  increment  of  data 

sec 

3 

1-10 

V 

FIO.O 

Velocity 

ft/sec 

11-20 

DIA 

FIO.O 

Diameter 

in 

21-30 

AX 

FIO.O 

Axial  moment 

Ib^-in 

31-40 

RHO 

FIO.O 

Density 

siug//r 

41-  50 

DEL 

FIO.O 

Canted  angle 

deg 

4 

1-70 

TIME 

9F8.3 

Time  table  for  spin. 

(Only  when  DT=  0) 

sec 

The  number  of  entry 

must  equal  to  NT. 

5 

1-72 

PI 

9F8.3 

Spin  table.  The  number 

of  spin  entry  must  equal 

rev/sec 

to  NT. 

6 

1-10 

AC(1) 

FIO.O 

Guessed  initial  spin 

11-20 

AC(2) 

FIO.O 

Guessed  C, 

> 

21-30 

AC(3) 

FIO.O 

Guessed  C, 
r 

31-40 

AC(4) 

FIO.O 

Guessed  C, 

ft 


<  N OTE>  Card  number  4  is  only  required  when  the  time  increment  (DT)  is  zero, 
otherwise  do  not  input  the  card. 

Card  number  6  is  only  an  optional.  If  you  don’t  have  guessed  value,  put  0.0  for 
each  entry. 
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APPENDIX  C 

PRCX3RAM  LISTING  AND  INPUT  EXAMPLE 
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PROGRAM  SPINPIT( INPUT, TAPE5-INPUT, OUTPUT, TAPE6-OUTPUT, 

+  PUNCH, TAPE7-PUNCH, TAPE!, TAPE2,TAPE3,TAPE4  ) 

C  EQUATION  SOLVED  IS  PDOT-  CN2  +  CN3*P  +  CN4*P*P 
C  FTl-CONSTANT  TERM  OF  BEARING  FRICTION  (FT- LB)  (NEGATIVE) 

C  FT2-LINEAR  TERM  OF  BEARING  FRICTION  (FT-LB/RPS)  (NEGATIVE) 

C  CLPZ  -(PD/2V)  CLP  FOR  BARE  BODY  WITHOUT  FINS  (NEGATIVE) 

C  DEL  -  FIN  CANT  ANGLE  IN  RADIANS 
C  NT  -  NUMBER  OF  POINTS 

C  DT  >  TIME  INCREMENT  FOR  SPIN  VALUES  INPUTTED 
C  AC(1)  -  INITIAL  SPIN  VALUE 

C  AC(2)  =  ESTIMATED  CLD 

C  AC(3)  -  ESTIMATED  CLP.  MAY  BE  LEFT  BLANK 

C  AC(4)  -  ESTIMATED  CLP2 .  MAY  BE  LEFT  BLANK 

C  P-SPIN  RATE  (INPUT  IN  RPM) 

C  FTR  =  FACTOR  TO  CORRECT  CN(J) 

C  VMACH  -  MACH  NO 

C  DIA  “  DIAMETER  IN  INCHES 

C  AX  -  AXIAL  INERTIA  -LBF-INSQ 

C  TO  -  TEMPERATURE  -  DEG  F 

C  PO  -  BAROMETRIC  PRESSURE  -  PSI 

COMMON  CN(S) ,P( 200) ,0(5,200) ,TIME( 200) 

COMMON  NT, DTI 

DIMENSION  PI(200) , R( 200 ) , WXX( 5) , C( 5, 6 ) , AC( 5 ) 

DIMENSION  DIDENT(8) 

1000  FORMAT(8A10) 

1001  FORMAT(I5,F5.0) 

1002  FORMAT (5F1 0.3) 

1003  FORMAT(9F8.3) 

1004  FORMAT(//25X,  34HCLP  TOTAL  -  CLPZERO  +  CLP2* ( PD/2V )//) 

1005  FORMAT(/20X,  26HCALCULATED  INITIAL  SPIN  -  ,F7.2,  19H  REV/SEC  ERR 
COR  -  ,F4.2) 

1006  FORMAT(/36X,  lOHCLPZERO  "  ,F8.3,10X,  8HERROR  -  ,P8.6) 

1007  FORMAT(/30X,  16HCLDELTA*DELTA  =■  ,F9.4,  9X,  8HERROR  -  ,F8.6) 

1008  FORMAT(/39X,  7HCLP2  -  ,F8.3,  lOX,  8HERROR  -  ,  F8.6) 

1009  FORMAT(/21X,25HSTEADY  STATE  SPIN  RATE  -  ,  F9.4) 

1011  FORMAT(4F10.3) 

2000  FORMAT(lHl) 

2001  PORMAT(//20X,I5,2F12.2) 

2002  PORMAT(//10X,8F12.4/) 

2003  FORMAT(/10X,10P10.3) 

2004  FORMAT(//10X,I2,2F20.8) 

2005  FORMAT (/I OX, 5E1 8. 6) 

2006  FORMAT(/10X,2F15.4) 

2007  FORMAT(/10X,7F15.4) 

20 1C  FORMAT ( //20X, 4 HTIME , 8X , 8HOBSERVED , 7X, 1 OHCALCULATED , 6X , 8HOBS-CALC , 
C5X,9HCLP  TOTAL) 

2011  FORMAT(20X,  5H(SEC) ,8X,36H(ROLL  RATE  -  REVOLUTIONS  PER  SECOND)) 

2012  FORMAT ( ///I 5X,12HTOTAL  POINTS, 3X,10HDELTA  TIME) 

2013  FORMAT(///15X,8HDIAMETER,6X,4HAREA,5X,10HAX.  MOMENT, 2X, 8HVELOCITY, 
C  2X,12HDYN. PRESSURE, 3X,8HFIN  CANT/ 

C  17X,4H(FT),6X,7H(SQ-FT),4X,10H(LB-SQ-FT),2X,8H(FT/SEC) ,3X, 

C  10H(LB/SQ-FT),3X,9H(RADIANS)) 

2014  FORMAT(///10X,39HSPIN  INPUT  DATA  -  REVOLUTION  PER  SECOND/) 

2015  F0RMAT(1H1/4X,13HITERATI0N  NO. , 6X, 9HRESIDUALS , 7X, 14HPROBABLE  ERROR 
C  /20X,12H(REV/SEC)**2,11X,9H(REV/SEC) 

C  //10X,42HCORRECTION  OF  CONSTANTS  FOR  EACH  ITERATION//) 

2016  FORMAT (///I OX, 17HINITIAL  GUESSES— /18X, 2HP0, 9X, 3HCLD, 9X, 3HCLP, 

C  9X,4HCLP2) 

2017  FORMAT(///10X,25HTIME  INPUT  DATA  -  SECONDS/) 

2100  FORMAT(lH0/10X,8A10) 

10  READ  (5,1000)  DIDENT 

IF  (BOF(5)  .NE.  0.0)  GO  TO  300 

READ(5,1001)  NT, DTI 

DT-DTI 

DO  40  1-1,4 

CN(I)-0.0 
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AC<I)-0.0 
WXX(I)-0.0 
DO  20  J-1,NT 
Q(I,J)-0.0 
20  CONTINUE 
40  CONTINUE 
DO  50  I-1,NT 
P(I)=0.0 
PI(I)=0.0 
R(I)=0.0 
50  CONTINUE 

READ(5.1002)  V,DIA,AX,RHO,DEL 
IF(DTI  .NE.  0.0)  GO  TO  60 
READ(5,1003)  (TIME(  I ) ,  I=>1  ,NT) 

60  READ(5,1003)  (PI( I ) , I-1,NT) 
READ(5.1011)  (AC(I), 1-1,4) 
QQ-0.5»RHO*V*»2 
Al-6. 2831852 
DEL-DEL*(Al/360.0) 

AIX-AX/144.0 

D-DIA/12.0 

S-3.14159*D*D/4.0 

WRITE(6,2000) 

WRITE(6,2100)  DIDENT 
WRITE(6,2012) 

WRITE(6,2001)  NT,DT 
WRITE(6,2013) 

WRITE(6,2002)  D, S, AIX,V,QQ,DEL 
IP(DTI  .NE.  0.0)  GO  TO  70 
WRITE(6,2017) 

WRITE(6,2003)  (TIME( I ) , I-l ,NT) 

70  WRITE(6,2014) 

WRITE(6,2003)  (PI( I) , I-1,NT) 
WRITE(6,2016) 

WRITE(6,2002)  (AC( I ) , I-l , 4 ) 
WRITE(6,2015) 

AXX-AIX/32.174 

CD-QQ*S*D/AXX 

DV-D/(2*V) 

CV-CD*DV 

DO  90  K-1,NT 

PI(K)-PI(K)*A1 

90  CONTINUE 

IP(AC(2)  .NE.  0.0)  AC(2)-AC(2)*DEL 
IF(AC(1)  .EQ.  0.0)  AC(1)=PI(1) 
IF(PI(NT-3)  .LT.  PI(1))  GO  TO  119 
IF(AC(2)  .NE.  0.0)  GO  TO  113 
IF(DT  .NE.  0.0)  GO  TO  112 
DO  91  1-1,4 

AA-(PI(I+1)-PI(1))/(TIME(I+1)-TIME(1)) 

CN(2)-CN(2)+AA 

91  CONTINUE 

AC(2)-(1.0/(4.0*CD) )»CN(2) 

GO  TO  113 

112  DO  92  1-1,4 

AA-(PI(I+1)-PI(1) )/(FLOAT(I)»DT) 
CN(2)-CN(2)+AA 

92  CONTINUE 

AC(2)-(1.0/(4.0*CD) )»CN(2) 

113  CONTINUE 

IF(AC(3)  .NE.  0.0)  GO  TO  114 

ICHECK-NT-5 

PRATIO-PI {NT)/PI ( 6 ) 

XXX-ALOG(PRATIO) 

TT-ICHECK*DT 

IP{DT  .EQ.  0.0)  TT-TIME(NT)-TIME(6) 
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AC(3)-XXX/(CV*TT) 

GO  TO  114 

119  CONTINUE 

IP(AC(3)  .NE.  0.0)  GO  TO  114 
ICHECK=NT-3 

PRATIO=PI(ICHECK)/PI(l) 

XXX=ALOG(PRATIO) 

TT-ICHECK*DT 

IF(DT  .EQ.  0.0)  TT=TIME(ICHECK)-TIME(1) 
AC(3)-XXX/(CV*TT) 

114  CONTINUE 

CN(1)=AC(1) 

CN(2)=AC{2)*CD 

CN(3)=AC(3)*CV 

CN(4)=AC(4)*CV»DV 

NR=0 

E=0.0 

100  SR=0.0 
CALL  INT 
DO  110  K=1,NT 
R(K)=PI(K)-P(K) 

SR=SR+R(K)**2 
110  CONTINUE 
NR=NR+1 
E1=E 

E=0. 6745»SQRT(SR/(NT-4) ) 

SRR»SR/A1*»2 

EE=E/A1 

WRITE(6,2004)  NR,SRR,EE 
IP(ABS(El-E)-lE-04)  200,200,120 

120  IP(NR-15)  130,130,200 
130  DO  140  I»l,4 

DO  140  J=l,5 
C(1,J)=0.0 
140  CONTINUE 

DO  150  1=1,4 
DO  150  J=l,4 
DO  150  K=1,NT 

C(I,J)=C(I,J)+Q(I,K)*Q{J,K) 

150  CONTINUE 

DO  160  1=1,4 

DO  160  K=1,NT 

C(I,5)=C(I,5)+R(K)»Q{I,K) 

160  CONTINUE 

CALL  INV(C,4,5,WXX) 

DO  170  M=l,4 
CN(M)=CN(M)+C(M,5) 

170  CONTINUE 

WRITE(6,2005)  (C( I , 5) , 1-1 , 4 ) 

GO  TO  100 
200  CONTINUE 

AC(1)»CN(1)/A1 

AC(2)=CN(2)/CD 

AC(3)=CN(3)/CV 

AC(4)=CN(4)/(CV*DV) 

DO  210  1-1,4 

210  WXX{I)-SQRT(ABS(C{I,I)))*E 
WXX(1)-WXX(1)/A1 
WXX(2)-WXX(2)/CD 
WXX(3)-WXX(3)/CV 
WXX(4)-WXX(4)/(CV»DV) 

C4C-SQRT(ABS(CN(3)**2-4*CN(2)»CN(4) ) ) 
PST-(-CN(3)-C4C)/(2*CN(4) )/Al 
PSS-(-CN(3)+C4C)/(2*CN(4))/Al 
PINT-PI  (t‘T)/Al 
PST1-ABS( PST-PINT) 
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PST2-ABS( PSS-PINT) 

IF(PST1  .GT.  PST2)  PST-PSS 
WRITE(6,2000) 

WRITE(6,2100)  DIDENT 
WRITE(6,1004) 

WRITE(6,1005)  AC(1),WXX(1) 

WRITE(6,1006)  AC(3),WXX(3) 

WRITE(6,1008)  AC(4),WXX{4) 

WRITE(6,1007)  AC<2).WXX(2) 

WRITE(6,1009)  PST 
C 

CLPP2=-(AC(2)/(DV*PST*A1) ) 

WRITE(6,767)  CLPP2 

767  FORMAT(//20X,30HTHE  TOTAL  SPIN  DAMPING  MOMENT  , 

C  38HCOEPFICIENT  AT  STEADY  STATE  SPIN  IS  =  ,FI0. 
C 

WRITE(6,2000) 

WRITE(6,2IOO)  DIDENT 
WRITE(6,2010) 

WRITE(6,2011) 

DO  220  1=1, NT 
T=DT»( I-l) 

IF(DTI.EQ.O.O)  T=TIME(I) 

DPA=DV»P( I) 

PI(I)=PI(I)/A1 
P(I)=P(I)/A1 
R(I)=R(I)/A1 
CLP=AC( 3 ) +AC( 4 ) *DPA 
TCLP=CLP*QQ*S*D»DPA 
WRITE(6,2007)  T,PI(I) ,P(I) ,R(I) ,CLP 
220  CONTINUE 
GO  TO  10 
300  STOP 
END 

SUBROUTINE  INT 

COMMON  CM( 5) ,P(200), 0(5,200), TIME(200) 

COMMON  NT, DTI 
DT-DTI/5 
N=1 
1-6 

Q(l,l)=1.0 

P(1)=CN(1) 

100  IP(DTI  .EQ.  0.0)  DT={TIME(N+1)-TIME(N) )/5 

101  AKI=DT»(CN(2)+CN(3) *P(N)+CN(4)»P(N) »»2) 
PMl=P(N)+AKl/2 

AK2-DT*(CN(2)+CN(3)*PM1+CN(4)*PM1**2) 

PM2=P(N)+AK2/2 

AK3=DT» (CN( 2) +CN(3 ) *PM2+CN( 4 ) »PM2**2 ) 

PM3=P(N)+AK3 

AK4=DT* ( CN ( 2 ) +CN ( 3 ) *PM3+CN ( 4 ) *PM3 •»  2 ) 

P(N+l)-P(N)+(AKl+2*(AK2+AK3)+AK4)/6 

CP1=CN(3)+2*CN(4)»P(N) 

CP2-CN ( 3 ) +2*CN ( 4 ) ‘PM! 

CP3=CN(3)i-2*CN(4)*PH2 

CP4=CN(3)+2*CN(4)*PM3 

BK1-DT»CP1*Q(1,N) 

BK2=DT»CP2*(Q( l,N)+BKl/2) 

BK3=DT*CP3*(Q( l,N)+BK2/2) 

BK4-DT»CP4*(Q(1,N)+BK3) 

Q(1,N+1 )-Q(l,N)+(BKl+2*(BK2+BK3}+BK4)/6 

BK1=DT*(CP1»Q(2,N)+1.0) 

BK2=DT*(CP2*(Q(2,N)+BKl/2)+1.0) 

BK3-DT»(CP3*{Q(2,N)+BK2/2)tl.O) 

BK4-DT*(CP4*(Q(2,N)+BK3)+1.0) 

Q(2,N+l)-Q(2,M)+(BKI+2*(BK2+BK3)+BK4)/6 

BK1-DT*(CP1*Q(3,N)+P(N)) 
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BK2-DT* ( CP2 ♦ ( Q ( 3 , N ) +BK1/2 ) +PM1 ) 

BK3-DT*(CP3*(Q(3,N)+BK2/2)+PM2) 

BK4=DT*{CP4*(Q<3,N)+BK3)+PM3) 

Q(3,N+l)=Q(3,N)+(BKl+2*(BK2+BK3)+BK4)/6 

BK1=DT«(CP1*Q(4,N)+P(N)»*2) 

BK2=DT*(CP2*(Q(4,N)+BKl/2)+PMl**2) 

BK3=DT* ( CP3  * ( Q ( 4 , N ) +BK2/2 ) +PM2  *  *2 ) 

BK4=DT*{CP4*(Q(4,N)+BK3)+PM3**2) 

Q(4,N+l)=Q(4,N)+(BKl+2»(BK2+BK3)+BK4)/6 

N=N+1 

IP{N  .LT.  I  )  GO  TO  101 

I»I+1 

N=N-4 

P(N)»P(N+4) 

Q(l,N)=Q(l,N+4) 

Q(2,N)=Q(2,N+4) 

Q(3,N)«Q(3,N+4) 

Q(4,N)=0(4,N+4) 

IP(N  .LT.  NT)  GO  TO  100 
200  RETURN 
END 

SUBROUTINE  INV(C,NC,NCS1/WXX) 

DIMENSION  C(5,6) ,WXX(5) ,PIVOT{2),CC(5,10) 
3000  PORMAT(/10X,  20HDET  IS  EQUAL  TO  ZERO) 
NCT=NC*2 
NCP1=NC+1 
DO  10  1=1, NC 
DO  10  J=1,NC 
10  CC(I,J)=C{I,J) 

DO  20  1=1, NC 
DO  20  J=NCP1,NCT 
20  CC(I,J)=0. 

DO  30  1=1, NC 
30  CC(I,NC+I)=1. 

DO  205  1=1, NC 
PIVOT(l)=CC(I,I) 

DO  200  K=1,NC 
PIVOT(2)=CC(K,I) 

126  IF(K-I)  135,130,140 
130  DO  150  J=1,NCT 

IF(PIVOT(l))  134,210,134 
134  CC(K,J)=CC(I,J)/PIVOT(l) 

150  CONTINUE 
GO  TO  200 

IF(PIVOT(l))  136,160,136 
136  CC(K,J)=CC(K,J)-CC(I,J)*PIVOT(2)/PIVOTCl ) 
160  CONTINUE 
GO  TO  200 

140  DO  170  J=1,NCT 

IF(PIVOT(2))  145,170,145 
145  CC(K,J)=CC(K,J)/PIVOT(2)-CC(I,J) 

170  CONTINUE 
200  CONTINUE 
205  CONTINUE 
GO  TO  250 
210  write  (6,3000) 

250  DO  300  1=1, NC 
DO  300  J=1,NC 
300  C(I,J)=CC(I,J+NC) 

NCS-NCSl 

350  IF(NCS-NC)500,500.400 
400  DO  420  1=1, NC 
WXX(I)-C(I,NCS) 

420  C(I,NCS)-0. 

DO  450  J-1,NC 

450  C(I,NCS)-C(I,NCS)+C(I,J)»MXX(J) 
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NCS-NCS-1 
GO  TO  350 
500  CONTINUE 
RETURN 


END 

♦EOR 


SIMM, 

M819  RP 

SMOKE  CARTRIDGE  (  MACH 

=  0.3  ) 

25 

0.0  0.0 

335.55 

3.189 

0.504  0.002377 

0.057 

4.19 

1.94 

2.41 

2.76  3.05  3.29 

3.51 

3.72 

3.92 

4.51 

4.82 

5.10  5.43  5.80 

6.15 

6.48 

6.78 

7.08 

7.40 

7.76 

15.0  21.5  25.5 

30.5 

32.5 

1.82 

2.55 

3.24  3.87  4.44 

4.62 

5.00 

5.00 

5.85 

6.49 

6.67 

7.27  7.83  8.37 

9.00 

9.47 

10.00 

10.59 

10.91 

11.43 

20.25  23.90  25.25 

26.10 

26.10 

0.0 

0.0 

0.0  0.0 
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APPENDIX  D 

PROGRAM  OUTPUT  EXAMPLES 
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81MM,  M819  RP  SMOKE  CARTRIDGE  (  MACH  -  0.3  ) 


TOTAL  POINTS  DELTA  TIME 

25  0.00 

DIAMETER  AREA  AX.  MOMENT  VELOCITV  DVN. PRESSURE  FIN  CANT 

(FT)  (SQ-PT)  (LB-SQ-FT)  (FT/SEC)  (LB/SQ-FT)  (RADIANS) 

.2658  .0555  .0035  335.5500  133.8177  .0010 


TIME  INPUT  DATA  -  SECONDS 


1.940 

2.410 

2.760 

3.050 

3.290 

3.510 

3.720 

3.920 

4.190 

4.510 

4.820 

5.100 

5.430 

5.800 

6.150 

6.480 

6.780 

7.080 

7.400 

7.760 

15.000 

21.500 

25.500 

30.500 

32.500 

SPIN 

INPUT  DATA  -  REVOLUTION  PER 

SECOND 

1.820 

2.550 

3.240 

3.870 

4.440 

4.620 

5.000 

5.000 

5.850 

6.490 

6.670 

7.270 

7.830 

8.370 

9.000 

9.470 

10.000 

10.590 

10.910 

11.430 

20.250 

23.900 

25.250 

26.100 

26.100 

INITIAL  GUESSES— 

PO 

CLD 

CLP 

CLP  2 

0.0000  0 

.0000 

0.0000 

0.0000 
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ITERATION  NO.  RESIDUALS  PROBABLE  ERROR 

(REV/SEC) **2  (REV/SEC) 

CORRECTION  OP  CONSTANTS  FOR  EACH  ITERATION 


1 

41219.66284045 

29.38302266 

.183458E+01 

-.966689E+00 

-.581463E-01 

.428149E-04 

2 

2858.58084748 

7.86950737 

-.314311E+01 

.328480E+01 

-.733116E-01 

.102561E-03 

3 

174.12605513 

1.94224578 

.325840E+00 

-.592481E+00 

.215574E-01 

-.209033E-03 

4 

14.68677214 

.56407322 

.105435E+01 

-.171272E+01 

.537446E-01 

-.328188E-03 

5 

.85974779 

.13647649 

.341713E+00 

-.631507E+00 

.201828E-01 

-.119405E-03 

6 

.33173285 

.08477476 

.246289E-01 

-.438376E-01 

.138279E-02 

-.801743E-05 

7 

.33021108 

.08458009 

.391805E-03 

-.346121E-03 

.746384E-05 

-.353107E-07 

8 

.33021106 

.08458008 
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81MM,  M819  RP  SMOKE  CARTRIDGE  (  MACH  -  0.3  ) 
CLP  TOTAL  ■=  CLPZERO  +  CLP2*(PD/2V) 


CALCULATED  INITIAL  SPIN  = 

1.89  REV/SEC 

ERROR  = 

.05 

CLPZERO  = 

.004 

ERROR  = 

.000801 

CLP2  = 

-.183 

ERROR  = 

.010131 

CLDELTA»DELTA  - 

.0006 

ERROR  - 

.000013 

STEADY  STATE  SPIN  RATE  « 

26.7534 

THE  TOTAL  SPIN  DAMPING  MOMENT  COEFFICIENT  AT  STEADY  STATE  SPIN  IS  - 


-.00866 
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SIMM,  M819  RP  SMOKE  CARTRIDGE  (  MACH  -  0.3  ) 


TIME 

(SEC) 

OBSERVED 

(ROLL 

CALCULATED 

RATE  -  REVOLUTIONS 

OBS-CALC  CLP  TOTAL 
PER  SECOND) 

1.9400 

1.8200 

1.8898 

-.0698 

.0026 

2.4100 

2.5500 

2.6902 

-.1402 

.0023 

2.7600 

3.2400 

3.2881 

-.0481 

.0020 

3.0500 

3.8700 

3.7842 

.0858 

.0018 

3.2900 

4.4400 

4.1949 

.2451 

.0016 

3.5100 

4.6200 

4.5711 

.0489 

.0014 

3.7200 

5.0000 

4.9299 

.0701 

.0013 

3.9200 

5.0000 

5.2711 

-.2711 

.0011 

4.1900 

5.8500 

5.7307 

.1193 

.0009 

4.5100 

6.4900 

6.2734 

.2166 

.0007 

4.8200 

6.6700 

6.7965 

-.1265 

.0004 

5.1000 

7.2700 

7.2664 

.0036 

.0002 

5.4300 

7.8300 

7.8163 

.0137 

-.0000 

5.8000 

8.3700 

8.4274 

-.0574 

-.0003 

6.1500 

9.0000 

8.9993 

.0007 

-.0006 

6.4800 

9.4700 

9.5324 

-.0624 

-.0008 

6.7800 

10.0000 

10.0114 

-.0114 

-.0010 

7.0800 

10.5900 

10.4846 

.1054 

-.0013 

7.4000 

10.9100 

10.9827 

-.0727 

-.0015 

7.7600 

11.4300 

11.5343 

-.1043 

-.0017 

15.0000 

20.2500 

20.1429 

.1071 

-.0057 

21.5000 

23.9000 

24.0043 

-.1043 

-.0074 

25.5000 

25.2500 

25.1993 

.0507 

-.0079 

30.5000 

26.1000 

26.0038 

.0962 

-.0083 

32.5000 

26.1000 

26.1951 

-.0951 

-.0084 
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